Maximally-localized Wannier Functions in Antiferromagnetic MnO within the 
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We have calculated the maximally-localized Wannier functions of MnO in its antiferromagnetic 
(AFM) rhombohedral unit cell, which contains two formula units. Electron Bloch functions are 
obtained with the linearized-augmented-plane-wave method within both the LSD and the LSD+t/ 
schemes. The thirteen uppermost occupied spin-up bands correspond in a pure ionic scheme to 
the five Mn 3d orbitals at the Mni (spin-up) site, and the four O 2s/2p orbitals at each of the Oi 
and O2 sites. Maximal localization identifies uniquely four Wannier functions for each O, which 
are trigonally-distorted sp 3 -like orbitals. They display a weak covalent bonding between O 2s/2p 
states and minority-spin d states of Mn2, which is absent in a fully ionic picture. This bonding is 
the fingerprint of the interaction responsible for the AFM ordering, and its strength depends on the 
one-electron scheme being used. The five Mn Wannier functions are centered on the Mni site, and 
are atomic orbitals modified by the crystal field. They are not uniquely defined by the criterion 
of maximal localization and we choose them as the linear combinations which diagonalize the r 2 
operator, so that they display the D3d symmetry of the Mni site. 

PACS numbers: 71.15.-m,71.15.Ap,75.30.Et,75.50.Ee 
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INTRODUCTION 



The mean-field one-particle description of the elec- 
tronic struture of periodic crystalline solids is usually 
based on extended Bloch functions (BFs). Within the 
Born-von Karman periodic boundary conditions, the 
cyclic translational subgroup, which commutes with the 
effective one-electron Hamiltonian, contains N transla- 
tions by corresponding direct lattice vectors R. This 
Abelian subgroup has TV one-dimensional irreducible 
representations, which are labeled with wavevectors k 
within the first Brillouin zone (BZ). Therefore, extended 
Bloch states ^mk(r) in this description are classified 
with two quantum numbers, the band index m and the 
crystal-momentum k, and are obtained by diagonaliza- 
tion of the effective one-electron Hamiltonian. An al- 
ternative description can be derived in terms of local- 
ized Wannier functional! (WFs) w n (r — R), which are 
defined in real space via a unitary transformation per- 
formed on the Bloch functions. They are also labeled 
with two quantum numbers, the orbital index n and the 
direct lattice vector R indicating the unit cell they belong 
to. In contrast to BFs, WFs are useful, for example, in 
visualizing chemical bonds or in describing the dielectric 
properties of non-metallic materials. They can be con- 
sidered as the generalization to solids of the concept of 
"localized molecular orbitals" for finite systemsa. How- 



ever, one major problem in practical calculations within 
this representation is the non-uniqueness of WFs, related 
to the phase arbitrariness of the BFs, and to the arbitrary 
unitary transformations which can be performed on the 
BFs at any given k-point in the BZ. As demonstrated re- 
cently by Marzari and Vanderbilttl, this non-uniqueness 
can be resolved in principle by imposing the further con- 
dition of maximal localization. However, some residual 
arbitrariness remains, related to the choice of the local- 
ization criterion. 

The late transition-metal (TM) monoxides MnO, FeO, 
CoO, and NiO challenge the theory of electronic states 
since several decades. These highly correlated materials 
feature a Mott insulator character, and the conventional 
local-spin-density (LSD) schemes gets into difficulties 
with the localized TM d-orbitals, predicting in partic- 
ular incorrectly their spectral weight and their energy 
relative to O 2p states. All these oxides have the rocksalt 
structure in their paramagnetic phase, while below the 
Neel temperatuce Tn, a type 11 antiferromagnetic (AFM) 
ordering occurscl. The experimentally well-documented 
compound MnO, which is considered to be in the in- 
termediate charge-transfer/Mott-Hubbard regime, is a 
particularly suitable case study, since the available the- 
oretical schemes apply best to this material. Indeed, 
because of the exchange stabilization of its half-filled 
d shell, even LSD predicts its insulating character, al- 
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though the corresponding energy gap and magnetic mo- 
ment are much smaller than experimental data. Several 
electronic structure calculations have been performed for 
MnO, using traditional band-stcucture schemes, like un- 
restricted Hartree-Fock (UHF)Q and local-spin density 
(LSD), as well as more innovative approaches taking into 
account at various levels the large value of the on-sit. 
Coulomb repulsion, for the metal 3d states, like the SICl 
and the LSD. i UB methods, and the more recent model 
GW schemeErtj. Among the quantities calculated within 
the latter approach, special attention has been devoted to 
the quasi-particle spectrumel, and recently to the zone- 
center optic phonon frequencies and the Born effective 
charge tensorQ. Comparison with experimental data of 
the various physical quantities computed in the different 
schemes, has demonstrated a monotonic trend with the 
separation energy between occupied and empty d bands, 
which is too small within LSD, too high within UHF, 
and has about the correct value within model GW. This 
energy separation is strongly related to the on-site in- 
teraction U. The analysis of all these results, has been 
performed in terms of the extended Bloch states. 

On the other hand, no ab-initio investigation exists of 
the electronic states in AFM MnO in terms of localized 
Wannier functions. Keeping in mind the inherent limita- 
tions of the mean-field one-particle approaches, it would 
be however instructive to materialize within the WFs 
description how superexchange (which is responsible for 
AFM ordering) manifests itself in a first-principles one- 
electron picture. Furthermore, thq-anaximal localization 



II. METHOD 



A. FLAPW calculations 



method of Marzari and VanderbiltE 3 
the past to several periodic systems 
a low (trigonal) symmetry solid like AFM MnO, nor to 
a compound with partially filled cation d shell. 



iias-bcen applied in 
OO, but not yet to 



Motivated by the above reasons, we present in this 
work an ab-initio calculation of the maximally-localized 
Wannier functions of AFM MnO corresponding to the 
uppermost occupied bands, using the all-electron, full- 
potentiaLJinearized-augmented-plane-wave (FLAPW) 
methodEjiij. Because the trends in the electronic proper- 
ties from UHF to LSD are already known for this system, 
most of the computations have been performed for con- 
venience within the latter scheme. To investigate the ef- 
fects related to the on-site Coulomb interaction, we have 
also used the LSD + U method, implemented with the 
same FLAPW technical ingredientsEj used in the LSD 
computations. 

The manuscript is organized as follows. In Sec. II 
we give technical details of the FLAPW implementation, 
and structural informations relevant to AFM MnO. The 
necessary ingredients of our WFs calculations are then 
introduced. As we follow closely the formalism of Rcf. [|, 
we will not describe in detail the method itself. A general 
result regarding maximal localization of s-p WFs is de- 
ferred to an Appendix. In Sec. Ill, we present and discuss 



the results of our maximally-localized WFs calculations. 



Finally, in Sec. IV we draw our conclusions. 



All quantities presented and discussed in this paper 
have been_computed using the scmirclativiatic FLAPW 
methocOEj, expanded with local orbitalsEEl (lo) where 
appropriate. Inclusion of lo's in addition to the nor- 
mal FLAPW basis enforces mutual state orthogonality 
and increases variational freedom. This allows to treat 
the semicore Mn 3s, 3p states together with the va- 
lence states, and helps in dealing with the linearization 
of Mn 3d and O 2s, 2p states. Core states are cal- 
culated fully relativistically and self-consistently in the 
crystal potential. For the LSD calculations, the Hcdin 
and Lundqvist exchange-correlation functional has been 
adopted. In the calculations performed within the LSD 
+ U method, the values of the Hubbard and exchange 
constants, U — 6.9 eV and J = 0.86 eV, have been 
taken from Ref. [|. The atomic-sphere radii for Mn and 
O are chosen to be 2.0 and 1.8 a.u., respectively, and the 
FLAPW basis size is set to include all plane waves with 
energy up to 16.0 Ry. Four special points inside the ir- 
reducible wedge of the BZ were used for evaluating the 
charge density during self-consistency cycles. 



B. Structural details 

AFM ordering occurs in MnO below Tn = 117 K along 
any one of the (111) directions. The transition is accom- 
panied by a small crystallographic distortion that_txans- 
forms the cubic structure into a rhombohedral oneO: the 
90° angle between the lattice vectors increases only by 
0.624°. We have showrfl earlier that the effect of this 
distortion on the calculated zone-center optic phonon 
frequencies and Born effective charge tensor is negligi- 
ble compared to the anisotropy induced solely by the 
magnetic order. Therefore, throughout the present study 
we use the ionic positions of the perfect rocksalt chem- 
ical cell, with the experimental lattice constant value 
a = 4.435 A. In the AFM phase, the magnetic mo- 
ments of MnO are arranged into ferromagnetic sheets 
which are parallel to (111), while the direction of mag- 
netization in neighboring planes is reversed. The mag- 
netic and crystalline structure of MnO is shown in Fig. YU. 
The rhombohedral magnetic cell, whose volume is twice 
that of the paramagnetic rocksalt one, corresponds to the 
space group Dg d (R3m). Several choices for the prim- 
itive translation vectors are possible. However, in or- 
der to minimize discrete mesh effects in k-space integra- 
tions, and also for taking advantage of symmetry prop- 
erties (see below), we consider here the following primi- 
tive translation vectors (in the cubic coordinate system): 
t x = a (1, §, 5), t 2 = a (i, 1, |), and t 3 = a (|, ±, 1). 
The two equivalent anions Oi and O2, whose site sym- 
metry is C3V, are located at positions ±a(^, 3, 5), while 
the two non-equivalent cations Mni (spin-up) and Mn2, 



3 



whose site symmetry is D3d, are at (0, 0, 0) and a(l, 1, 1). 
With our choice of origin on Mni , the two Oxygen sites 
are therefore related by spatial inversion. 



Mn 2 




Mn-| 



FIG. 1: Antiferromagnetic structure of MnO. Magnetic mo- 
ments are arrayed in ferromagnetic sheets (shaded areas) par- 
allel to (111) planes. Mni, Mn2, and Oi,2 sites are given by 
large black, large grey, and small open circles, respectively. 
The trigonal coordinate axes are also indicated with long- 
dashed lines. 



C. Maximally localized Wannier functions 

We consider the general case of a composite group of 
energy bands labeled by their band index m = 1, 2, . . . , /, 
and connected among themselves by degeneracies along 
high symmetry lines, but isolated from all other bands at 
lower and higher energy. The phase arbitrariness (gauge 
dependence) of each Vmk) with respect to the set of al- 
lowed k values propagates to the corresponding Wannier 
functions |w„r), labeled by n = 1,2,..., / and lattice 
vectors R. This resulting non-uniqueness manifests it- 
self trough an arbitrary unitary matrix Uml appearing in 
the most general transformation between BFs and WFs, 
which is given by 

M^f dke-* k - R Y, U ™n\^), (1) 

I J JBZ m 

where V is the real-space primitive cell volume. A direct 
consequence of this arbitrariness is that linear combina- 
tions of the WFs in the set {|w ra R.)} also form a suitable 
basis. 

The strategy of Ref . || is to pick out from this arbitrary 
choice of WFs, the particular set which is maximally 
localized according to some criterion. Once this crite- 
rion has been chosen, and the composite group of bands 



specified, the search for the set of maximally localized 
WFs becomes a problem of functional minimization in 

(k) 

the space of the matrices Umn- The selected functional, 
which measures the sum of the quadratic spreads of the 
WFs probability density, is given by 

n = J2[(r 2 )n-(r)li (2) 

n 

where (r)„ = (w n0 \r\w n0 ) = 

i V/(2n) 3 J BZ dk (u„ )k | V k |u n , k ) and (r 2 }„ = 

(w n o\r 2 \w n o) n . In these expressions, u n ^ (r) = 

e"' k ' r J2m Umn "0mk (r) is a periodic function which can 
be obtained from the Bloch functions of the composite 
group of bands. Practically—one calculates the BFs on 
equispaced Monkhorst-Packt3 meshes of k points in the 
unit cell BZ (whose volume is equal to the conventional 
BZ one), built on the reciprocal lattice vectors. The 
grids have been offset in order to include T. With our 
chosen cell geometry of MnO, the reciprocal lattice 
vectors are: gi = it/a (3, —1, —1), g2 = n/a (— 1, 3, —1), 
and g3 = 7r/a (— 1, —1, 3). We then express the matrix 
elements of the gradient Vk appearing in the localization 
functional in terms of finite differences. As shown in 
Rcf. [j| the only information needed is the overlap matrix 
M^', ' — (u n k|w n ' i k+b)j where b are vectors connecting 
each mesh point to its nearest neighbors. From these 
quantities, the gradient of the spread functional f2 with 
respect to an infinitesimal unitary transformation SUmn 
of the BFs can also be evaluated. Once this gradient 
is computed, the minimization can take place via a 
steepest descent or conjugate-gradient algorithm. 

It is important to stress that if the minimum of £1 is 
flat (as it will be shown to be the case for some of the 
MnO WFs), its precise location may be hindered by the 
various numerical approximations involved in the calcu- 
lations. In this work, we have especially taken care of the 
following features, which clearly become irrelevant in the 
limit of a dense mesh (N — > oo, b — > 0): (i) the regular 
mesh {k} has been chosen in order to have the lattice 
symmetry, i. e. it transforms into itself by application of 
the point group operations R of the crystal; (ii) the shell 
of b vectors used in the finite-difference formula for Vk 
has been constructed on the 12 vectors of type {332} and 
{1 I 0} (in internal units). They are the sides midpoints 
of cubes which can be built on particular points of the 
rhombohedral mesh, and this set of 12 points is therefore 
compatible with the highest possible symmetry. 

Finally, and following the procedure described in Eqs. 
(62)-(64) of Ref. [| we start the minimization procedure 
by constructing a set of "trial functions" in the unit cell, 
which are an initial guess of the final WFs. We use 
Gaussians centered on atomic sites and modulated by 
an appropriate combination of spherical harmonics, with 
a rms width value such that the gaussian is negligible 
outside the corresponding atomic sphere. A unitary ro- 
tation among the initial BFs is then performed in order 
to maximize their projection on the trial functions. 
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It has been shownt3 that Wannier functions can be cho- 
sen with well-defined symmetry properties, and forming 
(generally reducible) representations of the crystal point 
group. In order to investigate the nature and possible re- 
ductions of these representations in the case of MnO, and 
also to be in the position to use group-theory methods 
for further analysis, we need to calculate the represen- 
tation matrix elements (uvo|-Pr l^no)) where Pr is the 
transformation operator corresponding to the operation 
R of the point group, and the WFs are in the central cell. 
For simplicity, we restrict the formalism to symmorphic 



spacegroups (which is actually the case for D| d ). Using 
Eq. (|lj) with the k points discretization and R = 0, we 
have for the rotated Wannier function 

PR Ko) = ^EE U rnl P R IVW), (3) 

kesz m 

and, with the "periodic gauge" condition V>m _Rk (r) = 

"0m i?k+G k (r) 



(w n ' \P R \w n0 )NV - E E U rn^J U ml (VWk' |-Pfl|^mk)iW 

k,k' m.m' 

= jf E_ E +Gk) (VW Kk +Gk |PR|^k)v. (4) 

kEBZ™.™' 



The last braket is trivial only in the case of non- 
degenerate BFs. In the general case, on the other hand, 
we have 



Pr ipmk (r) = ipmk (R 1 r) 

= ^ D m , m (R, k) ip m > R u (r), 



(5) 



where D m i m (R 1 k) is a unitary transformation associated 
with the symmetry operation R, and the summation is 
over states degenerate with Vw flk (r) . If this latter state 
is non-degenerate, D m i m (R, k) = 5 m ' m , and Eq. (||) be- 
comes the usual formula. We note that the above formu- 
las allow to decompose \w n o) into a linear combination 
of basis functions transforming according to the relevant 
irreducible representations of the point group. 



III. RESULTS 

For our Wannier-function description of the electron 
states in AFM MnO, we have considered as a single com- 
posite group the topmost thirteen valence bands (O 2s, 
O 2p, Mn 3d). In this way, the corresponding Wannier 
functions display the highest localization, and therefore 
can be considered as the "elementary building blocks" of 
the occupied electron states. As it is convenient to in- 
terpret some of our results in terms of combinations of 
atomic orbitals, we give in Table || the relevant Oh repre- 
sentations corresponding to s, p, and d orbitals and their 
reduction into irreducible representations of and C3 V , 
together with their corresponding basis. Cartesian coor- 
dinates in the parent cubic reference system are written 
using x, y and z symbols, while X, Y and Z correspond 
to the trigonal reference system (Z axis along the cubic 



[111] direction). Throughout this work, we use the Bethe 
notation for the irreducible representationsE2l 



TABLE I: Classification of s, p, and d atomic orbitals accord- 
ing to the irreducible representations of Oh together with their 
reduction into irreducible representations-of D3d and C3 V , and 
the corresponding basis. Bethe notationEj is used. The cubic 
and trigonal coordinate systems are denoted by a;, y, z and 
X, Y, Z, respectively. 



Oh irreducible 
representations 



r 4 " 



r+ 



Reduction into (C3 V ) irreducible 
representations and basis 



IT (Tij 



r 2 " (ri) 



r 3 " (r 3 



r 3 + (T 3 ) 



rt (Ti 
T 3 + (Ta 



{Px + Py + Pz)/V3 = pz 
(p« +Py ~ 2p z )/v / 6 =px 

(Px ~Py)/V2 =PY 



(d yz + d zx + d xy )/\/3 = d Z 2 
(d yz + d zx - 2d x „)/v / 6 
(dyz - d zx )/V2 



We have first studied the convergence of the spread 
functional £1 in terms of the mesh density, using a sam- 
pling of the Brillouin zone BZ defined by the Monkhorst- 
Pack meshes v y. v x v with v = 2, 4, 6, and 8. Analysis 
has been done on the various terms of the localization 
functional O discussed in Ref. ||. We found that a satis- 
factory level of convergence is achieved using the 6x6x6 
mesh, which represents a good compromise between ac- 
curacy and computer burden. Also, such a mesh is large 
enough to prevent difficulties related to finite grids and 
occurence of periodic WFs replica. The 6x6x6 sampling 
has been retained throughout this work. 
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The minimization of the localization functional, which 
determines the unitary matrix U m n, has been performed 
using a mixed strategy: a simple fixed-step steepest- 
descent procedure is used during the first iterations, fol- 
lowed by several iterations with a conjugate-gradient 
procedure which is reset to steepest descent every 100 
iterations. Because of the delicate convergence of the 
Mn WFs, and in order to be on the safe side, we used a 
very large number (40000—80000) of iterations. 




FIG. 2: LSD Wannier functions of AFM MnO corresponding 
to the whole spin-up valence-band complex (13 bands), (a) 
The 2s/2p WF of Oi with center along the [111] direction; 
(b) one of the 3 equivalent Oi WFs, with center close to the 
[100] axis. The Mni, Mn2, and Oi,2 sites are denoted by large 
black, large grey, and small white spheres, respectively. Light 
grey and dark grey colours indicate the positive and negative 
amplitudes of the WFs contours. 

Several choices of trial orbitals have been tested to ini- 
tialize efficiently the minimization procedure. In particu- 
lar (see Section |n]) , we have considered the following sets 
of initial trial orbitals: (i) Gaussians times real combi- 
nations of spherical harmonics s, p x , p y , and p z on the 
Oi and O2 Oxygen sites, and d z 2, d x 2_ y 2, d xy , d yz , and 
d zx on the Manganese Mni site. These real harmonics, 
reported in Table H, form the basis of reducible represen- 
tations of C3V (O site) and of D 3 d (Mni site), and are 
oriented along cubic axes; (ii) same centers as above, but 
with Gaussians modulated by spherical harmonics com- 
binations s, p x , Py, Pz, and d z ?, d X 2_ Y 2, d X y, d Y z, 
and dzx, z-e. basis functions of the irreducible repre- 
sentations of C3V and D3d; (hi) combinations of the two 
above possibilities. All these sets of initial trial orbitals 
eventually lead, at convergence, to the same WFs, which 
are real, in agreement with the discussion in Ref. |3| How- 
ever, the set (i) corresponds to a much faster convergence 
than the other choices. The reason is directly related to 
the shape and orientation of the final O 2s/2p WFs, and 
this point will be discussed further below. 

We give in Table || the converged Wannier centers and 
spreads corresponding to the four Oxygen WFs at the Oi 



site for the spin-up chanel. 



and r z are the compo- 



nents of the Wannier center position with respect to the 
corresponding atomic site. We give also the total spread 
of these four WFs, and the sum of the correspond- 
ing centers whose interpretation will be discussed below. 
WFs related to the O2 site have the same ft and opposite 
r. The spreads corresponding to the five LSD WFs (9- 



TABLE II: Center and spread of the four Wannier functions 
for the Oi site corresponding to the uppermost LSD spin-up 
valence-bands of MnO. Center components r x ,r y , and r z are 
given with respect to the atomic site in the cubic reference 
system. Sum of WFs centers and the total spread are also 
displayed. Values in parentheses correspond to the LSD + U 
scheme. Considering the C3V symmetry of the Oi site, WF 1 
belongs to Fi while WFs 2, 3, and 4 are equivalent to each 
other. WFs 5-8 for the O2 site have the same fi and opposite 
r. 



WF 


r x (A) 


r v (A) 


rz (A) 


fi (A 2 ) 


1 


0.174 


0.174 


0.174 


0.751 




(0.186) 


(0.186) 


(0.186) 


(0.686) 


2 


-0.387 


0.036 


0.036 


0.773 




(-0.366) 


(0.049) 


(0.049) 


(0.704) 


3 


0.036 


-0.387 


0.036 


0.773 




(0.049) 


(-0.366) 


(0.049) 


(0.704) 


4 


0.036 


0.036 


-0.387 


0.773 




(0.049) 


(0.049) 


(-0.366) 


(0.704) 


Total 


-0.141 


-0.141 


-0.141 


3.069 




(-0.082) 


(-0.082) 


(-0.082) 


(2.797) 



13) at the Mni site are given in the first line of Table III. 
By symmetry, these WFs are centered on the Mni site 



TABLE III: The spreads (in A 2 ) of the five Wannier func- 
tions (9-13) centered at the origin (Mni site) corresponding 
to the uppermost LSD and LSD + U spin-up valence-bands 
of MnO. Total spreads are also indicated. The meaning of 
s-LSD and s-LSD + U is explained in the text. 





9 


10 


11 


12 


13 


Tot (9-13) 


LSD 


0.6220 


0.6218 


0.6497 


0.6493 


0.6493 


3.193 


s-LSD 


0.5950 


0.5950 


0.6248 


0.6899 


0.6899 


3.194 


s-LSD+[/ 


0.4713 


0.4713 


0.5133 


0.5163 


0.5163 


2.488 



(origin). Because the two Oxygen sites are equivalent, 
corresponding results for the spin-down chanel can be 
simply obtained by reversing Mni and Mn2. The total 
spread (sum of the 13 individual spreads) is 9.332 A 2 . 
We note also that the sum of centers has all cartesian 
components equal to zero. Indeed, this quantity repre- 
sents the electronic polarization (modulo a lattice vector) 
corresponding to valence states, and has to be zero as the 
system is centrosymmetric. In order to analyze in detail 
Table |fi|, it is useful to construct the matrix elements 
of the representation of the symmetry operations of the 
point group on the basis of the 13 Wannier functions. 
This is done by using Eq. (^). First, we have verified by 
constructing the multiplication table that these 12 matri- 
ces form indeed a representation of D3d . For non-trivial 
operations of the point group D3d, they display the fol- 
lowing structure: (i) there is an 8x8 block corresponding 
to O 2s /2p Wannier functions. Each line of this block 
consists of zeros, except for one element which is equal 
to 1. Because the origin is on the Mni site, the operations 
of D3d are simply transforming all WFs on both O sites 
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into each other; (ii) there is a 5x5 block corresponding to 
Mni 3d Wannier functions. This block does not display a 
well-defined structure, and in general all its elements are 
non-zero; (iii) the two off-diagonal blocks (connecting 
O 2s /2p and Mni 3d WFs) have all their elements equal 
to zero, at an accuracy better than 10 -7 . Therefore, the 
13 x 13 representation is block-diagonal, and in particular 
the 5x5 Mni 3d block makes a (reducible) representation 
by itself. 

We discuss first the O 2sj2p Wannier functions. The 
two non-equivalent WFs corresponding to site Oi are 
shown in Fig. |2j, with a viewpoint which is about the 
same as in Fig. [H The equivalent WFs corresponding to 
site O2 can simply be obtained by inversion symmetry 
with respect to the origin. In the absence of discrim- 
ination between Mn! and Mn 2 sites, the Oxygen site 
would have cubic (Oh) symmetry, and the 4 converged 
WFs originating from s and p states on the same Oi site 
would be sp 3 hybrids. This is a general result, which 
is demonstrated in Appendix from localized atom- 
centered s—p orbitals, the maximal localization algorithm 
always favors sp 3 -like Wannier functions, provided that 
the symmetry is cubic or higher. This situation occurs, 
e.g., in a hypothetical ferromagnetic (FM) MnO crystal, 
which will be briefly described later. In AFM MnO, Mni 
and Mn 2 sites are non-equivalent, so that the symmetry 
of the Oxygen site is reduced to C3 V . However, remi- 
niscence of the sp 3 hybrids persists: disregarding for the 
moment the d admixture, we see in panel (a) of Fig. ||that 
the WF 1, with Ti symmetry, displays clearly the shape 
of a non-symmetric p z orbital, with its center shifted 
from the atomic site along the [111] direction (j) z -s mix- 
ing), as also indicated in Tabic ||. The three other WFs 
are equivalent to each other under the operations of C3V 
(they have the same spread) , and their centers are located 
much closer to the cubic axes than to the diagonals. This 
can be explained by considering the interaction with the 
3d orbitals of Mn 2 : all these WFs display a substantial 
bonding 3d contribution from the three neighboring un- 
occupied Mn 2 sites. As the first-neighbor shell of Mn 
sites is an octahedron centered on each O, the above in- 
teraction favors the orientation of the equivalent WFs 2, 
3, and 4 along cubic axes. We display in panel (b) of 
Fig. U the WF 2 whose center is close to the [100] axis. 
The sum of WFs centers given in Table |l| is a measure 
of the trigonal distortion of the sp 3 -like orbitals, since 
it would be vanishing in a perfect cubic environment. 
This configuration of the converged WFs justifies also 
our choice of trial orbitals: p x , p y , and p z orbitals are 
already good guesses of final WFs 2, 3, and 4. During 
the iterations, most of the effort is devoted to building 
up WF 1 (the p z -s hybrid oriented along [111]). We 
note also in Fig. ^ small features on the neighboring Mni 
sites. They correspond to a weak bonding contribution 
from unoccupied states on these sites. 

We consider next the five Mni 3d Wannier functions. 
They are exactly centered at the Mni site and their 
spreads do not display any well-defined symmetry. When 
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FIG. 3: Same as Fig. [2[ but for the five Mm 3d WFs within 
LSD, after diagonalization of {r 2 )„i„. The WFs (a) and (b), 
as well as (c) and (d), are partners of a representation of 
D3d, while (e) belongs to Ff. 



convergence of O 2s /2p Wannier functions is practically 
achieved, further iterations of the localization process 
have the only effect of modifying slightly (esentially by 
amounts comparable to the numerical noise) the indi- 
vidual spreads of the five Mni 3<i WFs, while their to- 
tal spread remains unchanged. These five WFs are not 
uniquely defined by the criterion of maximal localization. 
In fact, any set of WFs obtained by applying a unitary 
transformation to the original set, will in general have 
different individual spreads, but the same value of the 
total spread. Since it is advisable that the final WFs dis- 
play explicitely as much symmetry as possible, we exploit 
this non-uniqueness and perform a symmetry reduction 
of the 5-fold representation of the Mni 3d WFs, which 
decomposes into Tf + 2 according to D3d symmetry. 

We note that the Hermitian operator r 2 , which appears 
in the localization functional Eq. (||), has the full rota- 
tional symmetry, and therefore its matrix representation 
{r 2 ) n 'n = (w„'o|r 2 |w ri o) belongs to the representation 
of D3d . These matrix elements can be calculated through 
a simple generalization of Eq. (23) of Ref. ||, and are given 
by 



{wn'o\r 2 \w n0 ) = Wb [2 5 n > 



(k,b) 



k.b 

(k,b) . 



M. 



(6) 



We choose the five new WFs as the linear combinations of 
the original set which diagonalize the r 2 matrix, Eq. (|^). 
We get in this way an extra unitary transformation which 
allows us to update the matrices The resulting 

WFs, which display the expected (r^ + 2 T^ ) symme- 
try, are indicated with the header "s-LSD" in Table III, 
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and are shown in Fig. |[ We note that an arbitrariness 
(up to a two-by-two unitary transformation) remains in 
the definition of the partners in the two T% representa- 
tions. The figure shows that the WFs 9-13 are essentially 
atomic orbitals modified by the crystal field. 

The above picture depicting the whole set of WFs in 
AFM MnO has been confirmed by similar calculations 
performed for a hypothetical FM MnO crystal. In this 
case, the symmetry is cubic (space group Ojj), and the 
majority-spin chanel contains 18 valence states (O 2s, 
O 2p, and Mn 3d), while the minority-spin chanel has 8 
valence states (O 2s, O 2p). As in AFM MnO, the Mn 3d 
WFs are found to be essentially atomic orbitals. By con- 
trast, the O 2s/2p WFs in both chanels are here undis- 
torted sp 3 hybrids (see Appendix |A|). In the minority- 
spin chanel, both lobes of each hybrid display a substan- 
tial bonding 3d contribution from the three neighboring 
unoccupied Mn sites, while these features are absent in 
the majority-spin chanel (occupied Mn sites). The sp 3 - 
like Oxygen WFs, with center along the [111] direction, 
are shown in Fig. [l] for both spin chanels. 
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FIG. 4: LSD Wannier functions of FM MnO. The sp 3 -like O 
WF with center along the [111] direction for (a) majority, and 
(b) minority-spin. Notations are the same as in Fig. ^. 

The covalent interaction in AFM MnO between O 
and Mn2 3d orbitals, clearly singled out in Fig. || for 
the case of the spin-up chanel (we remind that corre- 
sponding results for spin-down are obtained by reversing 
Mni and Mn2), is the fingerprint of superexchange in a 
one-electron picture. Superexchange mechanism results 
in fact from indirect many-body interactions, involving 
metal atoms (here Mm and Mn2) with an intervening 
Oxygen. The shared covalency of nearest-neighbor pairs 
of magnetic ions leads to an antiferromagnetic alignment 
of their moments. As a matter of fact, supposing, e.g. a 
spin-up WF on the Mni site, the neighboring Mn WF 
has the choice between two spin orientations, correspond- 
ing to a FM or an AFM ordering. The FM choice will be 
energetically less favourable, due to the required orthog- 
onalization to the Mni WF, while in the AFM alignment, 
the spin part of the wavefunction accounts for this neces- 
sary orthogonalization. The above textbook explanation 
of AFM ordering can be discussed in light of our calcu- 
lations, showing for the first time a real-space represen- 
tation of the WFs for an AFM system. The maximally- 
localized WFs bring out the relevant physical informa- 
tion, namely the covalent interaction with empty Mn 3c? 



states. We also notice that the FM alignment, on the 
other hand, is favoured by kinetic energy. This can be 
seen in Figs. [| and |] from the larger number of covalent 
interactions (6 instead of 3 in AFM case) associated with 
the O sp 3 -like WFs, which leads eventually to a broader 
band-width. 

It is useful at this point to contrast this rather sim- 
ple one-electron picture of AFM MnO, obtained in real 
space using the localized Wannier-function description, 
with the more entangled picture based on the energy 
band-structure of the extended Bloch-function descrip- 
tion. For this purpose, we show in the left panel of Fig. || 
the spin-up band structure of MnO along the X — T — Z 
directions, calculated within the LSD scheme. Starting 
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FIG. 5: Energy bands of antiferromagnetic MnO calculated 
along the X — V — Z lines of the rhombohedral magnetic zone 
using the local-spin-density approximation (a), and the LSD 
+ U scheme (b) . Shaded areas correspond to the fundamental 
gap. Irreducible representations at F (Bethe notation) are 
also given on the left and right of the energy axis, for states of 
primary Mn and O origin, respectively. Bands corresponding 
to states with high probability in the Mn atomic spheres are 
drawn with dashed lines. 

from the lowest energies, we find first a and a 
states at k = 0. They arise from Ti O 2s orbitals on 
the two O sites (related by inversion symmetry) , forming 
symmetric (T^) and antisymmetric (r^~) combinations, 
respectively. The primary origin of the next composite 
group of 6 bands is O 2p triplets on the two O sites, which 
are split into a Ti singlet and a doublet by trigonal 
C3V site symmetry. They lead to one set of symmet- 
ric combinations (r^, T^"), and one set of antisymmetric 
combinations (r^, T^) between the two sites. Next, the 
occupied t^g orbitals of the Mni atom are split at T into 
a doublet plus a singlet (r^ , Tf), while the e g orbitals 
remain a doublet (r+). These Mni 3d orbitals give rise to 
the uppermost occupied composite group of five bands. 
Clearly, states of the Mni 3d and O 2p complexes which 
belong to the same (or Tjj") representation interact 
with each other, so that eventually both sets (Mni d), 
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rf(0 p) and T J (Mni d), T^(0 p) consist of bonding and 
antibonding partners of Mnx 3d/0 2p hybrids. A mirror 
picture (relative to Fermi energy) of the above descrip- 
tion applies to the conduction bands which originate from 
the Mn2 3d orbitals, if we exclude the lowest conduction 
band at T, which is free-electron-like, has mostly Mn 
4s /O 3s character, and is not relevant to the present dis- 
cussion. The above description would leave completely 
empty the spin-up 3d states of the Mn2 atom. However, 
a covalent interaction (clearly depicted in Fig. || within 
the Wannier-function description) between O 2p and the 
unoccupied Mn2 3d states takes place, and some amount 
of spin-up charge appears at the Mn2 sites, which would 
be totally empty in a completely ionic picture. For in- 
stance, the highest antibonding valence state at T has 
60% d-occupancy in the Mni sphere, 7% p-occupancy in 
each of the two O spheres, and 14% d-occupancy in the 
Mn2 sphere. 

With the above results in mind, it is instructive to in- 
vestigate WFs within other one-electron schemes, lead- 
ing to different degrees of wavefunction localization, co- 
valent interaction, and magnetic moments. Data corre- 
sponding to the LSD + U scheme are presented in Ta- 
bles [n] and m for the spin-up chanel. Values for the 
Mni 3d WFs are those obtained after diagonalization of 
(r 2 ) n 'n- All individual spreads, as well as the total spread 
(whose value is 8.083 A 2 ), are smaller than in the LSD 
case: the O 2s /2p and Mn 3d Wannier functions are more 
localized by 9% and 22%, respectively, and this is consis- 
tent with the narrower valence band width. This prop- 
erty appears also in Fig. [g, where we show the LSD + U 
result for the WF 2 derived from O 2s /2p states. This 
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FIG. 6: (a) The 2s/2p WF of Oi as in Fig. | (b), but cal- 
culated within the LSD + U scheme; (b) charge density dif- 
ference (LSD + U minus LSD) corresponding to the WF of 
panel (a). Light grey (a "+" sign has been added for clarity) 
and dark grey colours indicate positive and negative contour 
values, respectively. 

LSD + U orbital is similar to the corresponding one cal- 
culated within LSD, and displayed in Fig. ||b. The most 
noticeable difference appears in the weaker covalent inter- 
action between O and the Mn2 3d orbitals. This is clearly 
visible in the right panel of Fig. || where we display the 
d— like differential charge density (plsb+u — Plsd) asso- 
ciated with this particular Oxygen WF. The five Mn! 3d 
WFs, obtained within LSD + U, are very similar to the 
LSD ones, and are essentially atomic orbitals modified 



by the crystal field. 

The more localized character of all WFs within the 
LSD + U scheme, together with the weaker covalent in- 
teraction between O and Mn2 3d orbitals, is the finger- 
print in the Wannier-function description of the larger 
on-site interaction U of cation 3d states. It is worth to 
realize that this simple trend is more easily detectable 
in terms of Wannier functions than Bloch functions. To 
this end, we show in the right panel of Fig. || the LSD + 
U band structure of MnO. It differs in several respects 
from the LSD one. The energy gap (2.0 eV), though still 
smaller than the experimental one (3.8 — 4.2 eV), is closer 
to it than the LSD value (0.8 eV). The energy separation 
between occupied and empty d states, which is due to the 
on-site interaction U between cation 3d states, has also 
increased. Furthermore, the three Mni 3d bands corre- 
sponding to the lower and states (t2 g ) do not form 
any more a composite group of bands by themselves, but 
rather mix with the O 2p complex. This results in a nar- 
rowing of the valence band width (6.3 eV, against 7.2 eV 
in LSD). 

The distinct features of the WFs computed within 
LSD and LSD + U can be used for discussing the be- 
havior of the Born dynamical charge tensor Z* within 
these different many-body approximations. In fact, as 
Z* can be defined in terms of the relative displacement 
of WFs centersEil, it carries along a relevant informa- 
tion. In particular, a p—d interaction between valence 
and_conduction states in mixed ionic-covalent materi- 
alsE3 leads to anomalous Z*, while atomic-like occupied 
WFs correspond to a nominal value. Therefore in MnO, 
one expects smaller \Z*\ (i.e. closer to 2, its nominal 
value) within the LSD + U scheme, to which correspond 
more localized WFs, than within LSD. We performed 
the calculation of the [111] component of this tensor, us- 
ing the WF description, and obtained within LSD the 
value \Zn \ = 2.36, which is consistent with the one ob- 
tained previouslyu using the Berry phase approach. LSD 
+ U , on the other hand, gives as expected a smaller value 
\Zn\ = 2.26. These results relate the behavior of two 
important physical quantities (the superexchange energy 
and Z*) to the same information contained in the AFM 
MnO WFs, namely the magnitude of their amplitudes on 
unoccupied Mn2 sites in the different approximations. 



IV. CONCLUSION 

We have used the maximum localization criterion for 
constructing Wannier functions of AFM MnO within 
both LSD and LSD + U schemes. This is the first appli- 
cation of the method to a magnetic low-symmetry sys- 
tem, with a partially filled d band. The LSD + U ap- 
proach has been considered here, as it improves the LSD 
description of the ground-state properties of strongly 
correlated materials. The Wannier-function description, 
which can be viewed as a generalization to extended sys- 
tems of the localized molecular-orbitals description of 
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molecules, provides in the case of MnO a simple picture 
complementary to the Bloch-function one. We obtain 
two distinct groups of WFs: (i) the five Mn 3d WFs cen- 
tered on the Mni site which are atomic orbitals modified 
by the crystal field; (ii) the four trigonally distorted sp 3 
WFs associated with each of the Oi and O2 sites, which 
are not centered on these sites. The latter display a sub- 
stantial covalent bonding with 3c? states on Mn2 sites, 
which is consistent with the AFM ordering of this ma- 
terial. Concerning the Mn WFs, they are not uniquely 
defined by the criterion of maximal localization. We have 
exploited this feature to obtain symmetrized maximally- 
localized WFs. We impose this further condition by di- 
agonalization of the r 2 operator. 
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APPENDIX A: MAXIMALLY LOCALIZED 
WANNIER FUNCTIONS FROM ONE S-LIKE 
AND THREE P-LIKE ORBITALS FOR CUBIC 
SITE SYMMETRY. 

We consider the 4-dimcnsional Hilbert space H4 
spanned by 4 orthonormal localized basis functions \<f>i), 
\(p2), |<^3), and |<^>4), centered on a given site (taken here 
as the origin) with cubic point-group symmetry (T, Tj,, 
O, Td, or Oh). \<j>i) is the basis function of a representa- 
tion Tt (for T, O, and T d ) or T+ (for T h and O h ). \<j> 2 ), 
|(/)3), and |04) are the basis functions of a threefold rep- 
resentation T 4 (for T and O), T 5 (T d ), (T h ), or 
(Oh). More precisely, | ^2 ) = \4>x) 6 x-row of the three- 
fold representation, and Ife) = C^ 1 \4> x ), \4>a) — C^\(j) x ), 
where C3 is the rotation by 27r/3 around the [1 1 1] 
axis. We consider real \4>i) functions, as the maximally- 
localized WFs turn out to be real, apart from an arbitrary 
overall phase factom. Finally, we suppose that the \4>i) 
(i = 1, ... ,4) are sufficiently localized so that (<f>i \ r \<f>j) 
and ((j>i \ r 2 \<pj) exist and are finite. 

We demonstrate below that minimizing within the 
localization functional, Eq. (||), leads to sp 3 -like Wan- 
nier functions, as opposed to the original atomic-like 
\4>i) functions. The solution is unique and corresponds 
to one of the two tetrahedra defined by the cubic point- 
symmetry group if the latter is T or Td- For the higher 
symmetry groups Th, O, and Oh, maximal localization 
leads to an infinity of sp 3 -like solutions constructed from 
\4>i) and unitary transformed three-fold functions \4>2)i 
03), \4>i) according to any of the 3x3 orthogonal matri- 
ces (a,/3, 7) representing a rotation of the cartesian 
co-ordinates of a vector. 



As the operator r 2 belongs to the identity representa- 
tion, it is diagonal in the \4>i) basis, and we have 



i\r 2 \4>i) = (£l s il p Cl p Hp), 



(Al) 



where cubic symmetry has been taken into account. Sim- 
ilarly, writing r = xi + yj + zk. and using projection 
operators, we find after some algebra 



Ai Aj Ak 

Ai 5k 63 

Aj 6k 6i 

Ak 63 6i 



(A2) 



with A = (0 s |£|(/>£), where £ = x, y, z, and <5 = 
{<t>x\y\4>z) = {4>y\ z \4>x) = (4>z\x\4> y ). The phase arbitrari- 
ness of the orbitals allows us to set A > 0. We note that 
symmetry requires i5 = in the case of Th, O, or Oh- The 
total spread of the original atomic-like \cj>i) functions is 

^4> = E^l r2 \4>i) ~ |(&| r \&)\ 2 ] = H s + 3H P . (A3) 



We consider now the four orbitals obtained from the 
\4>j) by a general orthogonal transformation O, given ex- 
plicitely by 



\lPi) = E ^ \</>j) (h 3 = 1) 



■ 4). 



(A4) 



There are 6 independent matrix elements, and 10 or- 
thono rma lit y co nditions for columns or rows. Using 
Eqs. (Al), (A4) and the orthonormality conditions, we 
obtain 

I>*i r2 = Et °ii n * + (°« + °* + °^ Q p i 

i i 
i 

= n s +3fi p , (A5) 



which verifies the trace invariance property of Eq. (Al) 
under orthogonal transformations, and gives for the total 
spread Cl^, of the functions defined by (|A4|) the inequality 



Hi! 



E 



|^)-^|(^|r|^ 

i i 

n s + m p -J2\{^\r\^)\ 2 



(A6) 



Therefore, the atomic-like \cj>i) functions correspond to 
the maximum value of the spread, and the problem of 
the total spread minimization reduces to minimizing the 
second term — J^. \ r \ipi)\ 2 of the loc aliz ation func- 
tional, which can be evaluated using Eqs. ( A2) and ( A4). 

We consider next the two particular orthogonal trans- 
formations O and O' of the \4>i) corresponding to the two 
possible choices for the four sp 3 hybrid orbitals and 
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\4>i) (written together below \ fpi^)), 
by 



and given explicitely where i = 1, 



i(|^)±|^)±|^ 3 >±|04» 
i(|0!> ± |&) T |0 3 ) T 104)) 
|(|&)=F|&>±|03>T|04» 



I ^m)-^( I ^)T | 02)T | 03) + | 04)) 



(A7) 



(111) 



(111). Using the O trans- 



formation and Eq. ( A2), we get for the \ r matrix 
elements 



-<5)(i+j-f 
(A - S) i 
(A-5)j 
(A - 5) k 



k) 



(A 



(A - 5) i 
f 5)(i-j-k) 
-(X-5)k 
-(A-5)j 



(A- 



(A 

-(A 

-(A-S)i 



<5)k 



(A - 5) k 
-(A-S)j 
-(A-5)i 

(A + «)(-i-j+k)y 



(A8) 



We note that taking the O' transformation corresponds 



to changing A into — A in Eq. ( A8) for the (^|r|^) matrix 
elements. Using this latter equation, together with (A£), 
we obtain the total spreads of the \ipt) and \"ip[) WFs 

% = E[^k 2 i^)-i(^i r i^)i 2 ] 



= fl s + 3Q P - 3(5 + X) 2 

= - 3(5 + x) 2 < n^, 



and 



3(5- A) 2 < fL 



(A9) 



(A10) 



Two cases have to be considered according as 5 is zero 
(point groups Th, O, Oh) or not (point groups T and 
T d ). 

If 5 = 0, the two spreads f2~ and fi~, are both equal 
to — 3 A 2 . We demonstrate first that this latter value 



is the mini mum of t he spread, Eq. (A6). To this end, 
using Eqs. ( |A2| ), ( A4), and the orthonormality conditions 

T o 2 



1 and 2 x = 1, we are lead to minimizing 

^ - r> -4A 2 £ oKo^ + ol + ol) 

i 

= f! -4A 2 + 4A 2 ^O^ (An) 



with the constraint 

= ^01-1 = 0. 



(A12) 



The conditions for an extremum of 0^,, introducing the 
Lagrangian multiplier A are 

°"" + A^- = 16A 2 Ol + 2A O n = 0, (A13) 



dOa dO a 
from which we get 



O 2 — Ci 1 — Ci 1 — ri 2 — — 
u n — u 2i — u 3i — u 4i — A i 



(A14) 



and the result 

minf^ = - 3A 2 = Q~ = Q~ < fi^. (A15) 

We demonstrate next that in addition to O and O', 
there is in fact an infinity of orthogonal transforma- 
tions producing orbitals with this same minimum 
spread value. For this purpose, we introduce the operator 
R (a, (3, 7) corresponding to a geometrical rotation speci- 
fied by its Euler angles, and its three-dimensional unitary 
representation (a, (3, rfmm 1 , whose basis functions 
are the spherical harmonics Y\ m , with m = 1, 0, — 1. 
Consistently with our choice of the basis functions |</>2), 
03), and |04), we consider instead the equivalent, real 
orthogonal representation (a,/3,7)^' with £, £' = 
x, y, z, which is the usual transformation of the carte- 
sian co-ordinates of a vector. We then build from this 
matrix the following orthogonal transformation acting in 
the 4-dimensional Hilbert space TLa 



7) = 



1 













£>« (a,)8,7)«' 








(A16) 



We write |^) = J2k ^ T,j V ( a > 0>l)kj \<f>j), and get 
after some algebra, using Eqs. (A2) and (A4) 



f2^> — 3 A^ 



Va, /?, 7. 



(A17) 



It is important to remember that the \<pi) are basis func- 
tions of irreducible representations of a cubic group, 
and not of the full rot ation group. Therefore, apply- 
ing the transformation (A16) on the original basis func- 
tions |0i) does not result in rotated 0i (r) functions, i.e., 
0i(i? _1 (a,/?,7)r), but in linear combinations of the orig- 
inal basis functions, defining the same Hilbert space H.4, 
and carrying the same chemical information as the orig- 
inal functions. 

When 5^0 (case of T and Tj), the two transfor- 
mations O and O' produce and orbitals with 
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different total spreads 



n~<n~, 
n~>o~ 



if 8 > (A S > 0) 
if S < (A 5 < 0) 



(A18) 



We prove finally that fi ~ for 5 > (or fi~ , for 6 < 0) 
is a local minimum of the total spread. To this end, 
we consider an infinitesimal orthogonal transformation 
of the (or |$)) 



0(e) 



J + e A- 



1 



e 2 A 2 



0(e 3 ), 



(A19) 



where A is an antisymm etric al operator defi ned b y 6 real 
parameters. Using Eq. (A8) and applying ( A19 ) to the 
sets of sp 3 hybrid orbitals \ipi) and |^), we obtain for the 
infinitesimal variations of the spreads, up to second or- 
der terms in e, positive semi-definite quadratic forms in 
the parameters of the A operator. Hence the correction 
to fi~ (or fi~ ) induced by the infinitesimal transforma- 

tion, are non-negative. The spread f2~ (£i^,) therefore 

corresponds for 6 > (5 < 0) to a local minimum of the 
localization functional. 



* Electronic address: postma@dpmail.epfl.ch 

1 G. H. Wannier, Phys. Rev. 52, 191 (1937). 

2 S. F. Boys, in Quantum Theory of Atoms, Molecules, and 
the Solid State, edited by P. O. Lowdin (Academic Press, 
New York, 1966), p. 253. 

3 N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 
(1997). 

4 K. Terakura, T. Oguchi, A. R. Williams, and J. Kiibler, 
Phys. Rev. B 30, 4734 (1984). 

5 W. L. Roth, Phys. Rev. 110, 1333 (1958). 

6 S. Massidda, M. Posternak, A. Baldereschi, and R. Resta, 
Phys. Rev. Lett. 82, 430 (1999). 

7 A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65, 1148 
(1990); Z. Szotek, W. M. Temmerman, and H. Winter, 
Phys. Rev. B 47, 4029 (1993). 

8 V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. 
B 44, 943 (1991). 

9 S. Massidda, A. Continenza, M. Posternak, and A. 
Baldereschi, Phys. Rev. Lett. 74, 2323 (1995). 

10 S. Massidda, A. Continenza, M. Posternak, and A. 
Baldereschi, Phys. Rev. B 55, 13494 (1997). 

11 N. Marzari and D. Vanderbilt, in First-Principles Calcu- 
lations for Ferroelectrics: Fifth Williamsburg Workshop, 
edited by R.E. Cohen (AIP, Woodbury, New York, 1998), 
p. 146. 

12 I. Souza, R. M. Martin, N. Marzari, X. Zhao, and D. Van- 
derbilt, Phys. Rev. B 62, 15505 (2000). 

13 H.J.F. Jansen and A.J. Freeman, Phys. Rev. B 30, 561 
(1984). 

14 S. Massidda, M. Posternak, and A. Baldereschi, Phys. Rev. 
B 48, 5058 (1993). 

15 A. B. Shick, A. I. Liechtenstein, and W. E. Pickett, Phys. 
Rev. B 60, 10763 (1999). 

16 D. Singh, Phys. Rev. B 43, 6388 (1991). 

17 J. B. Goodenough and A. Hamnett, in Zahlenwerte und 
Funktionen aus Naturwissenschaften und Technik, edited 
by O. Madelung, Landolt-Bornstein, New Series, Group 
III, Vol. 17g (Springer, New York, 1984), p. 201. 

18 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 
(1976). 

19 E. Kriiger, Phys. Status Solidi B 52, 215 (1972); 52, 519 
(1972). 

20 G. F. Koster, J. O. Dimmock, R. G. Wheeler, and H. Statz, 
Properties of the Thirty-two Point Groups (M.I.T. Press, 
Cambridge, 1963). 

21 R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 



1651 (1993). 

22 M. Posternak, A. Baldereschi, H. Krakauer, and R. Resta, 
Phys. Rev. B 55, R15983 (1997). 

23 A generalization of the method of Ref. ^ for the case of 
entangled energy bands has been given by I. Souza, N. 
Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 
(2001). 



